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ABSTRACT 


Military  electronics  are  frequently  operated  in  excessively  confined  spaces  aboard 
ships  and  aircraft.  This  limited  space  impacts  significantly  on  the  space  available 
for  cooling  equipment.  The  optimal  solution  is  the  development  of  one  universal, 
modular  rack  for  shipboard  or  aviation  use.  With  a  modular  design,  upgrades  to 
equipment  could  also  be  accompanied  by  an  upgrade  to  the  cooling  rack  itself  with 
very  little  additional  cost  or  difficulty.  A  water  cooled  avionics  rack  can  provide  suf¬ 
ficient  cooling  for  any  piece  or  combination  of  avionics  equipment  if  enough  water 
flow  paths  are  used,  the  water  is  at  the  appropriated  temperature  and  the  water 
is  properly  distributed  within  the  passages.  To  determine  if  the  cooling  medium, 
water,  is  sufficiently  distributed  within  a  modular  cooling  rack,  an  analysis  of  the 
flow  and  pressure  distribution  of  the  coolant  is  required.  This  thesis  presents  a 
computer  code  that  has  been  developed  as  an  initial  step  in  the  total  design  of  a 
modular  cooling  rack  for  avionics  equipment.  In  itself,  the  code  details  a  specific 
design  technique  and  allows  for  the  determination  of  whether  a  proposed  configu¬ 
ration,  including  source  location,  characteristics  of  the  cooling  water,  and  the  size 
and  shape  of  the  proposed  flow  passages  will  indeed  provided  a  proper  distribution 
of  the  coolant. 
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I.  INTRODUCTION 


A.  GENERAL 

This  thesis  describes  the  development  of  a  computer  code  to  analyze  the  flow 
and  pressure  distribution  of  water  in  an  avionics  cooling  rack.  This  code  allows  for 
variable  rack  design  and  flow  sources.  The  code  is  written  in  FORTRAN  77.  In 
addition,  the  code  can  be  run  on  any  IBM  or  IBM  compatible  personal  computer. 

B.  BACKGROUND 

Extremely  high  temperatures  are  a  primary  cause  of  avionics  equipment  un¬ 
reliability.  The  origin  of  the  thermal  problem  is  in  the  continuous  effort  to  develop 
lighter  and  smaller  components.  These  smaller,  more  densely  packed  components, 
by  virtue  of  the  heat  dissipation  within  a  small  volume,  frequently  operate  at  exces¬ 
sively  high  temperatures.  These  high  temperatures  result  in  component  derating, 
performance  degradation  and  accelerated  failure.  [Ref  1] 

Successful  thermal  management  of  electronic  systems,  under  development  in 
the  latter  part  of  the  1990s  will  require  the  removal  of  as  much  as  500  W  from  a 
single  chip,  at  heat  fluxes  between  50  to  100  W/cm^.  Volumetric  heat  release  rates 
can  also  be  expected  to  increase  dramatically  and  are  likely  to  exceed  10  W/cm^. 
Silicon  chips,  with  embedded  bipolar  circuits,  have  traditionally  been  maintained  at 
temperatures  ranging  from  65®  C,  for  commercial  computers,  to  between  110°  C  and 
125°  C  for  military  equipment.  [Ref  2] 

Due  to  necessity,  military  electronics  are  frequently  operated  in  confined 
spaces  aboard  ships  and  aircraft.  This  limited  space  impacts  significantly  on  the 
space  available  for  extensive  cooling  equipment.  Air-cooled  avionics  systems  are 
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currently  the  most  frequently  encountered  cooling  systems  used  in  military  appli¬ 
cations.  However,  it  is  common  practice  with  military  avionics  to  upgrade  the 
operational  capability  of  a  ship  or  aircraft  by  upgrading  only  one  or  two  pieces  of 
avionics  equipment  at  a  time.  With  an  air  cooled  system,  this  is  very  ineffective. 
By  replacing  equipment  in  one  of  these  systems  with  equipment  of  a  different  size 
and  shape,  the  cooling  airflow  of  the  original  design  is  disturbed  and  the  resulting 
airflow  around  the  new  (as  well  as  the  old)  equipment  is  less  than  optimal.  The 
original  cooling  system  was  designed  for  equipment  configured  in  a  certain  manner 
and  little  attention  is  paid  to  the  reoptimization  of  the  cooling  system  when  changes 
are  made.  Moreover,  cost  restrictions  are  also  an  important  factor  in  the  design  of 
military  cooling  systems. 

One  solution  for  the  military  applications  problem  is  the  development  of  a 
universal  rack  for  shipboard  or  aviation  use.  If  this  rack  were  modular,  it  would  allow 
for  a  great  amount  of  flexibility  in  design  and  provide  significant  cost  savings.  With 
a  modular  design,  upgrades  to  equipment  could  also  be  accompanied  by  an  upgrade 
to  the  cooling  rack  itself  with  very  little  additional  cost  or  difficulty.  By  developing 
one  universal  rack  for  all  military  applications,  a  significant  improvement  in  cooling 
systems  for  updated  designs  can  also  be  realized.  This  could  be  accompanied  by 
additional  cost  savings  through  the  elimination  of  unique  cooling  systems  for  every 
different  avionics  suite. 

The  rack  would  be  a  structure  that  could  accomodate  the  placement  of  several 
modular  electronic  packages.  One  might  imagine  a  “sfioeior”  configuration  with,  say, 
16  ^ holes”  in  which  16  electronic  packages  could  be  inserted.  Each  of  the  packages  is 
presumed  to  contain  electronic  components  mounted  in  a  variety  of  ways  (i.e  printed 
circuit  boards  or  on  brackets  attached  to  the  walls  of  the  package).  Such  a  structure 
would  possess  at  least  64  flow  paissages  exclusive  of  flow  passages  that  are  part  of 


the  fluid  source  or  pump.  The  objective  is  to  assure  that  each  flow  passage  carries 
enough  fluid  to  absorb  the  dissipated  heat  that  is  somehow  conducted  to  the  rack 
structure.  Of  course,  a  heat  exchanger  may  be  required  to  transfer  the  heat  to  the 
ultimate  heat  sink  which  may  be  the  environmental  air  or  sea  water. 

Thus,  the  basic  assumption  is  made  that  a  modular,  water  cooled  avionics 
rack  can  provide  sufficient  cooling  for  any  piece  or  combination  of  avionics  equipment 
if  enough  water  flow  paths  are  used,  the  water  is  at  the  appropriate  temperature  and 
the  water  is  properly  distributed  within  the  passages.  To  determine  if  the  cooling 
medium,  water,  is  being  sufficiently  distributed  within  a  modular  cooling  rack,  an 
analysis  of  the  flow  and  pressure  distribution  of  the  water  is  required.  This  can  be 
accomplished  using  a  computer  code  utilizing  node  analysis.  Such  a  computer  code, 
allowing  for  a  variable  number  of  branches  and  junctions,  is  presented  here  as  a  first 
step  in  the  development  of  a  universal  cooling  rack  for  military  applications. 

C.  BASIC  THEORY  AND  APPROACH 

1.  Node  Analysis 

The  analysis  of  the  cooling  rack  is  based  on  the  stipulation  that  any 
size  (diameter  and  length)  of  passage  may  be  used  in  the  construction  of  the  rack. 
Variables  in  the  program  include  the  density  and  viscosity  of  the  water,  the  number 
and  location  of  flow  sources,  the  number  of  water  paths  entering  each  junction,  the 
shape  of  the  flow  passages  (circular  or  rectangul  r)  and  the  length  and  equivalent 
diameter  of  each  passage.  The  purpose  of  this  degree  of  flexibility  is  to  allow  for  easy 
redesign  of  a  rack  in  the  event  that  it  does  not  meet  the  requirements  which  are  the 
proper  distribution  of  cooling  water  in  each  section  of  pipe.  By  varying  either  the 
rack  configuration  or  the  state  of  the  water  via  computer  input,  a  rack  that  provides 
the  proper  flow  distribution  can  eventually  be  proposed. 
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The  calculations  for  the  flow  in  the  passages  employs  a  matrix  oriented 
procedure  used  in  network  analysis.  The  network  analysis  approach  can  be  tailored 
to  flow  m  passages  by  proposing  an  analogy  of  the  "current”  to  fluid  flow  and  the 
"voltage”  to  the  fluid  pressure.  The  "resistance”  is  then  attributed  to  the  friction  in 
each  length  of  pipe.  Therefore,  each  length  of  pipe  will  have  a  resistanc'  associated 
with  it,  and  possibly  a  pressure  source  as  well,  depending  on  the  design  of  the  cooling 
rack. 

2.  Laminar  Flow 

The  computer  code  is  designed  to  calculate  a  laminar  flow  distribution. 
If  the  flow  is  in  transition  or  turbulent,  there  is  a  signiflcant  increase  in  the  amount 
of  frictional  resistance.  For  Reynolds  numbers  (Re)  less  than  or  equal  to  2100,  the 
flow  is  considered  laminar.  For  Reynolds  number  between  2100  and  10000,  the  flow 
is  in  transition  and  the  flow  is  turbulent  if  the  Reynolds  number  exceeds  10000.  The 
computer  code  calculates  the  Reynolds  number  for  each  flow  passage  and  provides 
a  warning  if  the  Reynolds  number  exceeds  2100. 

D.  SCOPE 

•  Chapter  II  explains  and  details  the  b2isic  code  required  to  analyze  the  flow 
and  pressure  distribution  of  water  in  a  proposed  cooling  rack  design. 

•  Chapter  III  describes  the  computer  code,  its  essential  capabilities  and  limita¬ 
tions,  associated  subroutines,  input  requirements  and  final  output. 

•  Chapter  IV  presents  the  results  from  several  case  runs  that  exhibit  the  flexi- 
bihty  and  capability  of  the  method. 

•  Chapter  V  concludes  with  future  development  efforts  and  the  application  po¬ 
tential  of  this  computer  code. 
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II.  AVIONICS  COOLING  RACK  PROBLEM 

FORMULATION 


A.  GENERAL 

The  flow  distribution  of  the  fluid  (in  all  branches)  of  the  cooling  rack  is 
determined  using  a  matrix  oriented  solution  technique.  The  general  equations  and 
the  matrix  solution  are  presented  in  what  follows. 

1  COOLING  RACK  BASIC  STRUCTURE 

The  basic  structure  of  the  cooling  rack  is  a  series  of  fluid  passages  fit¬ 
ted  together  using  standard  junction  components.  An  almost  unlimited  number  of 
passages  may  be  fitted  together.  An  example  of  one  possible  simple  combination  of 
passages  with  one  flow  source  is  presented  in  Figure  2.1. 


This  example  illustrates  how  variable  sized  flow  passages  may  be  used. 
The  rack  is  composed  of  b  straight  passages  called  branches  and  rif  junctions  called 
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nodes.  The  branches  are  numbered  because  the  computer  code  requires  that  basic 
information  for  each  flow  passage  branch  be  entered  separately.  The  nodes  are  also 
numbered  and  those  are  shown  in  circles.  Note  also  that  there  is  a  pressure  source 
located  at  branch-1.  For  convenience,  the  numbering  sequence  in  Figure  2.1,  begins 
with  a  1  on  the  upper  left  side  and  proceeds  to  the  right  and  then  down.  Although 
this  numbering  sequence  is  arbitrary,  maintaining  a  consistent  technique  facilitates 
later  adaptations  to  the  rack  and  comparison  with  other  rack  conflgurations. 

B.  VARIABLES  IN  DESIGN 

As  stated  in  a  previous  section,  the  length,  shape  and  the  diameters  of  the  flow 
p«issages  may  vary.  The  code  is  written  for  either  circular  or  rectangular  passages. 
For  rectangular  passages  an  effective  diameter  is  calculated.  Although  the  actual 
diameter  or  effective  diameter  may  vary  within  one  rack  design,  it  is  assumed  that 
either  all  circular  or  all  rectangular  passages  are  used  for  any  one  rack  configuration. 
All  length  and  diameter  size  information  is  entered  in  the  initial  part  of  the  program 
as  6  X  1  matrices  named  ELL(IB)  and  D(IB),  respectively.  The  density  and  the 
viscosity  of  the  cooling  water  are  input  as  the  variables  RHO  and  MU.  The  other 
input  and  output  variables  used  in  the  computer  code  are  summarized,  for  the 
reader’s  convenience,  in  Table  1. 


C.  DEVELOPMENT  AND  LINEARIZATION  OF 
FLOW  EQUATIONS 

1.  Laminar  Equations 

The  basic  equation  used  to  determine  the  change  in  pressure  or  pressure 
loss  within  a  fluid  passage  is  derived  from  the  D’Arcy  equation 


hf 


fLV^ 

29d, 


(2.1) 
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TABLE  2.1:  INPUT  AND  OUTPUT  VARIABLES  (FORTRAN) 


Variable 

Explanation 

A1 

Cross  sectional  height 

B 

Number  of  branches 

B1 

Cross  sectional  width 

D 

Branch  diameter  vector 

ELL 

Branch  length  vector 

IB 

Counter  for  branches,  IB  = 

MU 

Viscosity 

N 

Number  of  nodes 

NF 

Node  “from”  array 

NT 

Node  “to”  array 

PS 

Pressure  source  vector 

RHO 

Density 

In  this  equation,  hf  is  the  pressure  loss  due  to  friction  and  dg  is  the 
equivalent  diameter.  The  D’Arcy  equation  is  valid  for  steady  flow  within  passages 
running  full  of  liquid.  In  eq  (2.1),  /  is  a  friction  factor.  If  the  flow  is  laminar,  the 
friction  factor,  /,  as  shown  in  Figure  2.2,  can  be  represented  by  the  equation 


/  = 


Re 


(2.2) 


The  Reynolds  number,  by  definition  is 

Therefore,  /,  is  seen  to  have  an  alternate  form 

f  _  64m 
^  pVd, 

Substitution  of  eq  (2.4)  into  eq  (2.1)  yields 

.  32LVn 


(2.3) 


(2.4) 


(2.5) 
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Figure  2.2:  Moody’s  chart  for  the  friction  factor. 
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By  definition,  the  flow  rate,  Q,  is 


Q  =  VA  (2.6) 

The  area,  A,  of  a  circular  flow  passage  is,  of  course,  given  by 

A  =  -n^  (2.7a) 

4 

but  for  rectangular  passages  of  side  dimensions  a  and  6,  it  is 

A^ah  (2.7b) 

and  for  square  passages  where  a  =  6,  it  is 

A  =  a^  (2.7c) 

For  circular  passages,  d*,  is  simply  the  passage  diameter,  d,  and  for  rectangular 

passages,  dg  is  defined  as 


de  = 


CA 

P 


where  A  is  the  passage  flow  area  and  P,  for  channels  flowing  full,  is  the  passage 
wetted  perimeter.  In  order  to  make  the  equivalent  diameter  for  a  circular  passage 
equal  to  the  actual  diameter,  d,  C  =  4  so  that  for  a  rectangular  passage  with  side 
dimensions,  a  and  6,  the  equivalent  diameter  becomes 

2ab 


a  +  b 


(2.8a) 


In  the  event  that  the  passage  is  square  (a  =  6),  the  equivalent  diameter  becomes 


de  =  a 


(2.8b) 


Substitution  of  eqs  (2.7)  and  (2.8)  into  equation  (2.5)  yields  for  the 
circular  passage 


hf  = 


128LQfi 

gird*p 


(2.9a) 
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for  the  rectangular  passage 

8LtiQ{a  +  by 
^  /)<7(al>)3 

(2.96) 

and  for  the  square  passage 

Z2LfiQ 

(2.9c) 

pga* 

To  simplify  eq  (2.9a)  further,  let 

hj  =  P 

and 

R  - 

(2.10) 

gitd^p 

This  yields  the  equation 

P  =  RQ 

which  is  of  the  same  form  as 

V  =  RI 

Similar  adjustments  can  be  made  to  eqs  (2.9b)  and  {2.9c).  When  these  adjustments 

are  made,  it  is  seen  that  for  a 

rectangular  passage,  eq  (2.9b)  gives 

8LMa  +  by 

pg{aby 

(2.11) 

and  that  for  a  square  passage  eq  (2.9c)  provides 


^  _  32Lfi 
pga* 

Additional  losses  (other  than  passage  friction  already  discussed)  occur 
in  flow  passage  systems  and  cannot  be  disregarded  without  appreciable  error.  Com¬ 
pensation  for  entrance  and  exit  losses  must  be  considered  in  the  cooling  rack  system 
by  abiding  an  equivalent  length  of  straight  pipe.  The  equivalent  length  used  for 
square  edged  entries  is  20  diameters  (or  equivalent  diameters)  and  for  exits,  the 
equivalent  length  is  40  diameters  (or  equivalent  diameters)  [Ref.  2].  Therefore,  the 
additional  length  of  60  passage  diameters  accounts  for  both  exit  and  entrance  losses. 


10 


D.  NUMERICAL  SOLUTION  SCHEME 


We  first  consider  a  general  branch  (the  branch)  as  shown  in  Figure  2.3. 
The  branch  can  contain  a  pump  with  pressure  head,  with  an  associated  section 
of  pipe  of  resistance,  Rk-  The  branch  carries  a  flow,  qk^  and  possesses  a  total  head 
loss.  Pit.  Notice  that  the  orientation  of  the  head  loss  is  in  the  direction  of  the  flow 
(from  to 


Figure  2.3:  Branch  with  pressure  source. 


Also  observe  the  orientation  of  the  flow,  qk,  and  notice  that  by  the  analogy 
to  Ohm's  law  and  because  of  compatibility  (the  sum  of  the  losses  must  match  the 
head  available) 

Pk  =  Rk^k  +  Ptk 

This  may  be  written  for  all  k  branches  in  matrix  form 

P  =  RQ  +  P, 

Consider  a  pump  and  a  length  of  discharge  piping.  The  pump  can  be  repre¬ 
sented  as  a  flow  source  or  as  a  pressure  source  as  indicated  in  Figure  2.4. 
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Figure  2.4:  Alternative  source  arrangements  for  the  development  of  the 
Flow  source  4-4  Pressure  source  transformation. 


In  these  figures,  the  R's  represent  the  resistance  to  flow  of  the  discharge 
piping  and,  in  the  case  of  the  pressure  source,  the  +  and  —  indicate  the  discharge 
(+)  and  suction  (— ).  In  Figure  2.4a,  continuity  dictates  that. 


or 

p  =  Rgq  +  Rgqs 


In  Figure  2.4b,  compatability  (the  pressure  drops  around  the  simple  fluid  loop  must 
match) 

P  =  RvP* 


If  there  is  to  be  an  equivalence  between  the  two  (Figures  2.4a  and  2.4b)  then 


R  =  Rp  —  Rq 


and 


P*  =  Rq* 
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or  indeed 


R 

Thus,  when  we  represent  a  pump  as  a  pressure  source  with  a  resistance,  we  can 
immediately  transform  this  to  an  equivalent  flow  source.  The  computer  code  is 
written  for  the  source  input  to  be  a  flow  source. 

Note  also  that  continuity  dictates  that  for  the  branch 


Qk  =  ykPk  +  9sk 


And  for  all  k  branches  this  can  be  written  in  matrix  form  as 


Q  =  YP  +  Qs 


Next  consider  a  rack  containing  b  branches  and  nodes.  The  node  is  the 
datum  node  (the  node  at  the  suction  end  of  the  pump),  and  we  may  set  down  an 
Tit  X  b  augmented  node-branch  incidence  matrix  A“  with  elements 


{-f-1  if  branch  j  leaves  node  i 
—  1  if  branch  j  enters  node  i 
0  if  branch  j  is  not  incident 
upon  node  i 


(2.10) 


For  example  in  the  network  displayed  in  Figure  2.5  with  its  oriented  graph  shown 
in  Figure  2.6,  there  are  =  5  nodes  and  6  =  6  branches.  For  this  network,  the 
augmented  node  branch  incidence  matrix  A“  is 


A"  = 


0  0 

0  0 

1  0 

-1  1 

0  -1 


Note  that  for  each  column,  which  represents  a  single  branch,  there  should 
only  be  two  non-zero  entries,  a  1  and  a  —1.  These  entries  correspond  to  the  nodes 


Figure  2.5:  Network  for  example. 


Figure  2.0:  Oriented  graph  for  network  in  Figure  2.5. 


at  the  extremities  of  the  branch.  A  branch  that  leaves  a  node  shows  a  +1  while  a 
branch  that  leaves  a  node  shows  a  —1. 

Next,  we  form  a  node-branch  incidence  matrix  which  is  n  x  6.  This  is  done 
by  deleting  the  row  in  A“  that  corresponds  to  the  datum  node.  In  our  example, 

row-5  is  the  datum  node  and  thus 

■  1  1  0  0  0  0' 

0-1  1  1  00 

^  0  0-1010 
.0  0  0-1-11. 

The  vector  Q  is  defined  as  a  6  x  1  vector  representing  the  flow  in  all  k  branches 
and  its  elements  are  designated  as  The  product  of  A  with  Q,  AQ,  represents  a 
series  of  equations,  one  equation  for  each  node.  These  equations  represent  the  sum 
of  all  flow  in  and  out  of  the  respective  node.  Continuity  demands  that  product,  AQ 
be  null 

AQ  =  0 

and  this  simple  matrix  equation  is  a  statement  of  continuity  for  the  whole  system. 
To  illustrate 

'  <li  ' 

1  1  0  0  0  0  q2  {qi  "b  92) 

An—  0—1  1  1  0  0  93  _  (~92  +  93  +  94)  _  0 

^'^“00-1010  94  ~  (-93  +  95)  “  0 

0  0  0—1—11.  qs  .  (—94  —  95  +  96)  .  .  0  . 

.  96  . 

The  pressures  or  heads  at  the  nodes  are  represented  by  an  n  x  1  vector 
designated  by  H  and  the  branch  pressure  drops  are  represented  by  an  n  x  1  vector 
designated  by  P.  Reference  to  Figure  2.6  shows  that 

Pi  =  hi 
Pi  =  hi  —  Aj 
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P3  =  ^2  —  hs 

Pi  =  ^2  -  hi 

Ps  =  ^3  —  ^4 

P6  =  *4 

The  matrices  P  and  H  can  be  related  by  a  matrix  C 

1  0  0  0' 

1-1  0  0  [hr 

0  1-10  hj 

0  10-1  hs 

0  0  1-1  L  h4  . 

0  0  0  1. 

and  it  can  be  observed  that 

C  = 

Refer  now  to  the  general  branch  in  Figure  2.3  and  the  vector  representation 
for  the  flow  in  all  k  branches.  The  6x1  vector,  Q  which  represents  the  branch  flow 
rates  has  been  shown  to  be 

Q  =  YP  +  Q8  (2.11) 

Again,  we  carefully  note  the  orientation  of  the  flow  source.  Here  Y  is  an  n  x 
n  source  admittance  matrix.  Premultiply  eq  (2.11)  by  A  to  obtain 

AQ  =  AYP  +  AQs 

But,  P  =  A^H  and  AQ  =  0,  thus 

AYA'^H  +  AQs  =  0 

The  node  equations  are  then  formulated 

YnH  =  Q 


P  =  CH  = 
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where 


Yn  =  AYA^ 

and 

Q  =  -AQ» 

Therefore,  the  node  pressures  H,  can  now  be  determine  by 

H  =  Y-iQ 

Although  many  techniques  exist  for  the  solution  of  large  sets  of  linear  si¬ 
multaneous  algebraic  equations,  the  most  efficient  computationally,  appears  to  be 
the  Cholesky  reduction  followed  by  back  subsitution  that  is  employed  in  the  Gauss 
elimination  method.  The  only  restriction  on  the  use  of  Cholesky  reduction  is  that 
its  use  is  confined  to  symmetric,  positive  definite  matrices.  Here  it  is  fortunate  that 
Yn  is  always  positive  definite. 

The  Cholesky  reduction  is  based  on  the  premise  that  there  is  a  unique  lower 
triangular  matrix  that  permits  a  factorization  in  the  form  of  A  =  LL^  if  the  matrix 
A  is  symmetric  and  positive  definite.  Computationally,  the  Cholesky  reduction  is  a 
very  efficient  technique  in  that  it  only  requires  n{n  +  l)/2  storage  locations  for  L, 
rather  than  the  usual  n*  locations  required  by  other  methods.  [Ref.  3].  The  number 
of  operations  using  the  Cholesky  reduction  is  approximately  n^/6  rather  than  the 
usual  n^/3  required  for  most  other  decompositions  [Ref.  3]. 

The  Cholesky  method  to  solve  the  basic  system 


AX  =  B 


is  b2ised  on  finding  the  solution  to  two  equivalent  systems: 


L'^C  =  B  and  LX  =  C 
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The  elements  of  C  are  determined  by  the  algorithms 


Cl 


•Sll 


and 


bi  -  stiCi 


a  = 


/=i 


t  >  1 


Once  C  is  known,  X  can  be  found  using  back  substitution  as  employed  in  the  Gauss 


elimination  method,  [Ref.  5]  that  is 


Xn  = 


^nn 


and 


n 

Ci  -  ^  SuXi 
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III.  DESCRIPTION  OF  THE  COMPUTER 

CODE 

A.  PROGRAM  STRUCTURES  AND  CAPABILITIES 


1.  Restrictions  and  Limitations 

The  computer  system  used  is  an  IBM  or  IBM  compatible  personal  com¬ 
puter.  The  current  program  fixes  the  maximum  number  of  individual  fluid  passages 
allowed  at  300.  For  proper  use  of  this  program,  the  computer  must  have  a  minimum 
storage  requirement  of  2  M  Bytes.  A  detailed  computing  time  study  for  the  program 
has  not  been  undertaken  because  the  computing  time  changes  exponentially  with 
the  number  of  nodes  and  branches  selected  in  the  proposed  cooling  rack  design. 

The  computer  code  can  be  operated  using  the  following  information 

•  Laminar  Flow  Conditions 

•  Density  of  water  varying  from  x  to  y 

•  Viscosity  of  water  varying  from  x  to  y 

•  Up  to  100  branches 

•  Up  to  100  nodes 

•  Up  to  40  sources 

•  Any  diameter  (or  equivalent  diameter)  flow  passages 

•  Branches  of  any  length 

•  Metric  or  English  units 
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2.  Structure  of  the  Main  Computer  Code 

The  the  main  computer  code  including  the  three  associated  subroutines 
are  in  Appendix  A.  The  main  code  is  essentially  divided  into  two  major  sections, 
initialization/input  and  laminar  flow  solution  with  its  associated  output.  The  main 
function  of  the  initialization/input  section  are: 

1.  Set  up  the  problem  formulation  by  reading  in  the  necessary  values  (node  from, 
node  to,  length  and  cross  sectional  dimensions)  for  each  branch  of  the  cooling 
rack. 

2.  Read  in  the  viscosity  and  density  values. 

3.  Read  in  location  and  characteristics  of  each  pressure  source. 

The  second  section  calculates  the  pressure  and  flow  distribution  in  the 
rack.  It  includes  the  subroutines  MATMUL,  DECOMP  and  CHOLESKY.  The  main 
functions  of  this  section  are  to 

1.  Develop  the  node- branch  incidence  matrix,  A,  using  the  information  read  in 
section  one. 

2.  Calculate  the  effective  branch  lengths  to  account  for  losses  at  the  nodes. 

3.  Calculate  the  area  and  effective  diameter  of  each  flow  passage  cross  sectional 
area. 

4.  Calculate  the  branch  admittance  matrix,  Y. 

5.  Develop  the  flow  source  matrix,  Q#,  from  information  read  in  section  one. 

6.  Calculate  the  transpose  of  the  matrix  A. 
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7.  Calculate  the  node  admittance  matrix,  Yn- 


8.  Calculate  the  node  flow  source  vector,  Ig. 

9.  Calculate  the  node  pressure  vector,  H. 

10.  Calculate  the  inverse  of  Yn. 

11.  Calculate  the  branch  pressure  vector,  P. 

12.  Calculate  the  branch  flow  rate  vector,  Q. 

13.  Calculate  the  branch  Reynolds  numbers. 

14.  Provide  readouts  of  all  branch  flow  rates,  all  node  pressures  and  all  branch 
Reynolds  numbers. 

B.  DESCRIPTION  OF  SUBROUTINES 

1.  Subroutine  MATMUL 

This  subroutine  [Ref.  5]  multiplies  an  n  x  m  matrix  by  any  mx£  matrix 
to  form  an  n  X  /  matrix  and  is  called  whenever  matrix  multiplication  is  required. 

2.  Subroutine  DECOMP 

This  subroutine  [Ref.  6]  performs  the  decomposition  of  the  symmetric, 
positive  definite  matrix,  Yn  using  the  Cholesky  reduction.  It  is  called  by  subroutine 
CHOLESKY. 

3.  Subroutine  CHOLESKY 

This  subroutine  [Ref.  7]  provides  the  solution  of  a  linear  system  of 
equations  using  the  Cholesky  decomposition  method  for  positive  definite  matrices. 
This  subroutine  calls  subroutine  DECOMP. 
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IV.  RESULTS  AND  DISCUSSION  OF  CASE 

RUNS 


A.  CASE  RUN  ONE 

Case  one  was  run  using  the  configuration  shown  in  Figure  2.1  and  the  follow¬ 
ing  data: 

1.  Water  density  =  62.4  Ib/ft^ 

2.  Water  viscosity  =  8.8  x  10"“*  Ib/sec-ft 

3.  Basic  shoebox  design  with  8  nodes  and  12  branches 

4.  Length  of  end  branches  (1-4  and  9-12)  =  2.2  ft 

5.  Length  of  middle  branches  (5-8)  =  3.0  ft 

6.  Circular  passages  with  a  diameter  of  0.3  inches 

7.  One  flow  source  at  branch  one 

8.  Strength  of  source  was  0.4  gal/min 

The  results  from  case  one  are  in  Appendix  B.  It  should  be  noted  that  the 
Reynolds  number  in  each  branch  indicates  laminar  flow  and  that  there  is  cooling 
water  in  all  sections  of  the  sample  rack  design.  These  results  demonstrate  that 
the  analysis  technique  presented  here,  does  allow  for  the  determination  of  whether 
or  not  a  proposed  cooling  rack  design  will  indeed  provide  a  proper  distribution  of 
cooling  water.  It  also  verifies  that  the  laminar  flow  assumption  is  valid  for  some 
design  parameters. 


B.  CASE  RUN  TWO 

Case  run  two  demonstrates  the  flexibility  of  the  program  through  a  varia¬ 
tion  in  some  of  the  input  parameters.  The  basic  configuration  is  still  the  same 
configuration  shown  in  Figure  2.1. 

1.  Water  density  =  999.5  kg/m^ 

2.  Water  viscosity  =  13.1  x  lO”"*  kg/m-sec 

3.  Basic  shoebox  design  with  8  nodes  and  12  branches 

4.  Length  of  end  branches  (1-4  and  9-12)  =  1.0  m 

5.  Length  of  middle  branches  (5-8)  =  1.2  m 

6.  Rectangular  passages  with  a  height  of  1.3  cm  and  width  of  1.4  cm. 

7.  One  flow  source  at  branch  one 

8.  Strength  of  source  was  1.8  lit/min 

The  results  from  case  two  are  located  in  Appendix  B. 
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V.  CONCLUSION 


A.  GENERAL  COMMENTS 

This  computer  code  has  been  developed  as  an  initial  step  in  the  total  design 
of  a  modular  cooling  rack  for  avionics  equipment.  In  itself,  the  code  details  a 
specific  design  technique  and  allows  for  the  determination  of  whether  one  proposed 
configuration,  including  source  location,  characteristics  of  the  cooling  water,  and  the 
size  and  shape  of  the  proposed  flow  passages;  will  indeed  provide  a  proper  distribution 
of  the  cooling  water.  Without  proper  coolant  distribution,  the  cooling  rack  will  be 
inefficient  and  perhaps,  totally  ineffective. 

B.  ENHANCING  THE  MAIN  PROGRAM 
CAPABILITY 

The  next  step  in  the  development  of  a  computer  code  to  provide  more  flex¬ 
ibility  and  range  in  the  rack  design  is  the  inclusion  of  the  turbulent  flow  regime 
within  the  coolant  passages.  After  turbulent  flow  is  effectively  incorporated,  the 
final  phcise  in  the  total  rack  design  is  a  modification  for  heat  input  in  the  individual 
coolant  flow  passages. 
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APPENDIX  A 


PROGRAM  COOL 

c  Cooling  Rack  Calculations  Using  Node  Analysis 

r  INITIALIZATION 

r  MU:  fluid  viscosity  (Ib/ft-sec) 

r  RHO:  fluid  density  (Ib/ff^B) 

?  N:  number  of  junctions  in  system 

?  N1 :  number  of  junction  minus  the  datum  node 

?  B:  number  of  branches  in  system 

::  G:  acceleration  of  gravity  (ft/sec‘'2) 

r  PI :  pi 

INTEGER  N,N1,B, I, J,K,QB,S 
INTEGER  NT(40) ,NF(40) ,X, lUNITS 

REAL  MU,RHO,G,PI,A{40,40) ,D(40) ,ELL{40) ,R(40) ,MU1 
REAL  QS(40,1) ,AT(40,40) , QSl , II (40 , 1) , LENGTH 
REAL  E(40,1),P(40,1),Y(40,40}  ,YNI  (40,40) 

REAL  YP(40, 1)  ,Q(40, 1)  ,AY(40,40)  ,  YN(40,40) 

REAL  V(40, 1)  ,RE(40, 1)  ,al  (40)  ,bl  (40)  ,  AREA (40) 

REAL  IS (40, 1) , ISI (40, 1) ,RH01 
CHARACTER  ANS , ANSI , ANS2 

G=32.174 
PI=3 .1415926 
IUNITS=0 

^*«*******************************************«****«***************«****** 

c  THIS  PROGRAM  DETERMINES  THE  FLOW  RATE  AND  PRESSURE  DISTRIBUTION 

c  OF  COOLING  WATER  WITHIN  A  VARIABLE- SI ZED  AVIONICS  COOLING  RACK, 

c  VARIABLES  INCLUDE  RACK  DIMENSIONS  AND  WATER  CHARACTERISTICS. 

-^mt*********************************************************************** 


1040  FORMAT  (///'  ENTER  THE  TOTAL  NUMBER  OF  NODES  IN  THE  RACK  SYSTEM 
+  ' , 2X, \) 

1042  FORMAT  (BN,  13) 

1043  FORMAT  (//'  ENTER  THE  TOTAL  NUMBER  OF  BRANCHES  IN  THE  RACK  SYSTEM 
+  '.2X,\) 

222  WRITE (lOT, 1040) 

READ(IN,1042)  N 
WRITE (lOT, 1045)  N 

IF(N  .LT.  0  .OR.  N  .GT.  100)  THEN 
WRITEdOT,  6565) 

GO  TO  222 
END  IF 

223  WRITEdOT,  1043) 

READ (IN, 1042)  B 
WRITEdOT,  1046)  B 

IF(B  .LT.  0  .OR.  B  ,GT.  100)  THEN 
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WRITE (lOT, 6567) 

GO  TO  223 
END  IF 
c 
c 

1045  FORMAT  (//'  YOU  ENTERED  THE  NUMBER  OF  NODES  AS:',1X,I3\) 

1046  FORMAT  (//'  YOU  ENTERED  THE  NUMBER  OF  BRANCHES  AS ; ' , IX, 13 , / , \) 

6565  FORMAT  {/'  THE  NUMBER  OF  NODES  MUST  BE  GREATER  THAN  ZERO  AND  LESS', 

+  IX, THAN  100. ' ,\) 

6567  FORMAT  (/'  THE  NUMBER  OF  BRANCHES  MUST  BE  GREATER  THAN  ZERO  AND  LESS 

+  THAN  100 . ' , \) 
c 
c 

N1=N-1 
DO  5  1=1, B 
QS{I,1)=0. 

5  CONTINUE 

c 
c 
c 

55  CONTINUE 

WRITE (lOT, 7676) 

7676  FORMAT(/'  Do  you  want  to  work  in  British  units  or  SI  units  (B/S) ? 

+  '  f  \) 

READ (IN, 7677)  ANS 

7677  FORMAT (Al) 

IF(ANS  .EQ.  'B'  .OR.  ANS  .EQ.  'b')  GO  TO  7800 
IF (ANS  .EQ.  'S'  .OR.  ANS  .EQ.  's')  GO  TO  7678 
GO  TO  55 

7800  lUNITS  =  1 
c 

WRITE (lOT, 1000) 

1000  FORMAT  (//'  Input  Viscosity  (  x  10''4  Ibm/f t-sec)  '  , 2X, \) 

READ(IN,1003)  MU 

MU1=MU*. 0001 
WRITE (lOT, 1002) 

1002  FORMAT  (/'  Input  density  (lb/ft^3) ' , 2X, \) 

READ (IN, 1001)  RHO 

1001  FORMAT  (BN,  F8.4) 

1003  FORMAT  (BN,  F7.6) 

1004  FORMAT  (BN,  13) 

GO  TO  7801 

C 

C 

7678  WRITE (lOT, 7602) 

7602  FORMAT  (//'  Input  Viscosity  (Kg/m-sec  x  10*4)',2X,\) 

READdN,  1003)  MU 

WRITE (lOT, 7603) 

7603  FORMAT  (/'  Input  Density  (Kg/m^3) ' , 2X, \) 

READ (IN, 1001)  RHOl 

MUl-MU*. 00006719 
RHO=RH01*. 06243 
c 
c 

C  A:  NODE-BRANCH  INCIDENCE  MATRIX 
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C  INITIALIZE  A  MATRIX 
c 

7801  DO  20  1=1, N1 

DO  20  J=1,B 

20  A(I,J)=0. 

c 

DO  22  1=1, B 
NF{I)=0 
22  NT(I)=0 

c 
c 

c  DETERMININATION  OF  PASSAGE  SHAPE 
C 

c 

c 

8000  WRITE (lOT, 8001) 

8001  FORMAT (/'  IF  FLOW  PASSAGES  ARE  CIRCULAR,  ENTER  A  ZERO.',/, 

+  '  IF  PASSAGES  ARE  RECTANGULAR,  ENTER  A  ONE.',2X,\) 

READdN,  1042)  X 

IF(X  .EQ.  0)  THEN 
GO  TO  8004 
END  IF 

IF(X  .EQ.  1)  THEN 
GO  TO  8005 
END  IF 

8003  WRITE (lOT, 8006) 

8006  FORMAT (/'  ERRONEOUS  INPUT' \) 

GO  TO  8000 
C 
c 

c  INITIAL  DATA  INPUT  TO  DEVELOP  A, ELL  AND  D  MATRICES  FOR 
c  RECTANGULAR  PASSAGES 
c 

c  ELL :  LENGTH 

c  D :  DIAMETER 
C 
C 

8005  CALL  CLS 

IF  (lUNITS  .EQ.  0)  THEN 
DO  1493  J=1,B 
WRITE (lOT, 1051)  J 
WRITE (lOT, 1492) 

1492  FORMAT (/'  FROM  NODE,  TO  NODE,  LENGTH (m),  HEIGHT (cm),  WIDTH (cm)' 
+  /) 

READdN, 8202)  NF(J)  ,NT(J)  ,ELL(J)  ,al(J)  ,bl  (J) 

1493  CONTINUE 

7272  FORMAT (/'  YOU  ENTERED  THE  FOLLOWING  DATA.'\) 

7273  FORMAT (/'  BRANCH  FROM  NODE  TO  NODE  LENGTH (m)  HEIGHT (cm) ', 2X, 

+  'WIDTH (cm) ', /) 

7274  FORMAT (2X, 13 , 6X, 13 , 8X, 13 , 2X, F9 . 5, 4X, F7 . 5 , 4X, F7 . 5) 

7275  FORMAT(/'  IS  THIS  CORRECT  (Y/N) ?  '\) 

CALL  CLS 

16  WRITE (lOT, 7272) 

WRITE (lOT, 7273) 

DO  7277  J=1,B 

WRITE (lOT, 7274)  J,NF(J) ,NT(J) ,ELL(J) , al (J) ,bl (J) 
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7277  CONTINUE 

7276  WRITE (lOT, 7275) 

READ(IN,7677)  ANS 

IF{ANS  .EQ,  'Y'  .OR.  ANS  .EQ.'y')  GO  TO  7278 
IF{ANS  .EQ.  'N'  .OR.  ANS  .EQ.  'n')  GO  TO  7280 
GO  TO  7276 

7278  DO  7279  J=1,B 
ELL(J) =ELL(J) *3.28 
al(J)=al(J)*.3937 
bl(J)=bl(J)*.3937 

7279  CONTINUE 
GO  TO  1495 

7280  CONTINUE 

7281  FORMAT (/'  INPUT  ONE  BRANCH  NUMBER  TO  CHANGE  VALUES. '\) 

7284  WRITE (lOT, 7281) 

READ (IN, 1042)  K1 

7282  FORMAT(/'  INPUT  VALUES  FOR  BRANCH: IX, 13 , \) 

WRITE (lOT, 7282)  K1 

WRITE (lOT, 1492) 

READ (IN, 8202)NF{K1) ,NT(K1) ,ELL(K1) ,al (Kl) ,bl (Kl) 

7283  FORMAT (/'  DO  YOU  HAVE  ANYMORE  CHANGES  (Y/N)?'\) 

7285  WRITE (lOT, 7283) 

READ (IN, 7677)  ANSI 

IF(ANS1  .EQ.  'Y'  .OR.  ANSI  .EQ.  'y')  GO  TO  16 
IF (ANSI  .EQ.  'N'  .OR.  ANSI  .EQ.  'n')  GO  TO  7278 
GO  TO  7285 
END  IF 
C 

CONTINUE 

c 

c 

DO  8200  1=1, B 
WRITE (lOT, 1051)  I 
WRITE (lOT, 8201) 

8201  FORMAT(/'  FROM  NODE,  TO  NODE,  LENGTH(FT),  HEIGHT(IN),  WIDTH(IN) 

+  ',/) 

8202  FORMAT(2I3,F7.4,2F6.4) 

READ (IN, 8202)  NF(I) ,NT(I) ,ELL(I) ,al(I) ,bl(I) 

8200  CONTINUE 
c 

CALL  CLS 

15  WRITE (lOT, 7272) 

WRITE ( lOT  8207) 

8207  FORMAT (/'"  BRANCH  FROM  NODE  TO  NODE  LENGTH (FT)  HEIGHT ( IN) ', 2X, 
+  'WIDTH (IN) ', /) 

7374  FORMAT(2X,I3,7X,I3,7X,I3,4X,F9.5,6X,F7.5,3X,F7.5) 

DO  8277  J=1,B 

WRITE (lOT, 7374)  J,NF(J) ,NT(J) ,ELL(J) , al (J) ,bl (J) 

8277  CONTINUE 
8276  WRITE (lOT, 7275) 

READ (IN, 7677)  ANS 

IF (ANS  .EQ.  'Y'  .OR.  ANS  .EQ.'y')  GO  TO  1495 
IF(ANS  .EQ.  'N'  .OR.  ANS  .EQ.  'n')  GO  TO  8280 
GO  TO  8276 
8280  CONTINUE 
8284  WRITE (lOT,  7281) 
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READ (IN, 1042)  K1 
WRITE (lOT, 7282)  K1 
WRITE (lOT, 8201) 

READ (IN, 8202)NF(K1) ,NT(K1) ,ELL(K1) , al (Kl) ,bl (Kl) 

8285  WRITE (lOT, 7283) 

READ (IN, 7677)  ANSI 

IF(ANS1  .EQ.  'Y'  .OR.  ANSI  .EQ.  'y')  GO  TO  15 
IF (ANSI  .EQ.  'N'  .OR.  ANSI  .EQ.  'n')  GO  TO  1495 
GO  TO  8285 

1495  DO  1496  1=1, B 

IF  (NF(I)  .LT.  N)  THEN 
A(NF(I)  ,1)  =1. 

END  IF 

IF  (NT (I)  .LT.  N)  THEN 

A(NT(I)  ,1)  =  -1. 

END  IF 
al (I) =al (I) /12 
bl (I) =bl (I) /12 
AREA(I) =al (I) *bl (I) 

D(I)  =  (2*(al(I)*bl(I)))/(al(I)+bl(I) ) 

LENGTH=ELL ( I ) 

ADDL=60*d(I) 

ELL ( I ) =LENGTH  +  ADDL 

1496  CONTINUE 
GO  TO  8300 

C 

c 

C  RECEIVE  DATA  FROM  KEYBOARD  TO  DEVELOP  A,D  AND  ELL  MATRICES  FOR 
C  CIRCULAR  PASSAGES 
C 

1050  FORMAT (IX, 13 , 8X, 13 , 5X, F7 .4 , 4X, F6 . 4) 

9050  FORMAT (2 13, 2F7. 4) 

8004  CALL  CLS 

IFdUNITS  .EQ.  0)  THEN 

GO  TO  993 

ELSE 

DO  25  1=1, B 

1051  FORMAT  (/'  INPUT  THE  FOLLOWING  DATA,  SEPARATED  BY  COMMAS  FOR, IX, 
+  BRANCH: ' , I3,2X,\) 

1053  FORMAT  (/'  FROM  NODE, TO  NODE, LENGTH (FT) , DIAMETER (IN) ', 2X, /) 

WRITE (lOT, 1051)  I 
WRITE (lOT, 1053) 

READ (IN, 9050)  NF(I) ,NT(I) ,ELL(I) ,D(I) 

25  CONTINUE 

CALL  CLS 

18  WRITE (lOT, 7272) 

WRITE (lOT, 1059) 

1059  FORMAT (/'  BRANCH  FROM  NODE  TO  NODE  LENGTH (FT)  DIAMETER ( IN) ' , 

+  /) 

DO  9277  J=1,B 

WRITE (lOT, 1155) J,NF(J) ,NT(J) ,ELL(J) ,D(J) 

1155  FORMAT (2X, T3, 9X,I3,8X, 13 , 6X, F7 . 4 , 4X, F6 . 4 ) 

1255  FORMAT (2X, 13 , 7X, 13 , 7X, 13 , 5X, F7 . 4 , 5X, F6 . 4 ) 

9277  CONTINUE 
9276  WRITE (lOT, 7275) 

READ (IN, 7677)  ANS 


30 


9280 

9264 


9285 


c 

c 

993 

1055 


29 


19 

1077 


1277 

1276 


1280 

1284 


1285 


C 

111 


28 

c 
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IF(ANS  .EQ.  'Y'  .OR.  ANS  .EQ.'y')  GO  TO  26 

IF(ANS  .EQ.  'N'  .OR.  ANS  .EQ.  'n')  GO  TO  9280 

GO  TO  9276 

CONTINUE 

WRITE (lOT, 7281) 

READdN,  1042)  K1 
WRITE (lOT, 7282)  K1 
WRITE (lOT, 1053) 

READdN,  9050)NF(K1)  ,NT{K1)  ,ELL(K1)  ,D(K1) 

WRITE (lOT, 7283) 

READ (IN, 7677)  ANSI 

IF(ANS1  .EQ.  'Y'  .OR.  ANSI  .EQ.  'y')  GO  TO  18 
IF (ANSI  .EQ.  'N'  .OR.  ANSI  .EQ.  'n')  GO  TO  26 
GO  TO  9285 
END  IF 
GO  TO  26 


DO  29  J=1,B 
WRITE (lOT, 1051)  J 

FORMAT  (/'  FROM  NODE,  TO  NODE,  LENGTH (M) , DIAMETER (CM) ', /) 

WRITE (lOT, 1055) 

READdN,  9050)  NF(J)  ,NT(J)  ,ELL(J)  ,D(J) 

CONTINUE 

CALL  CLS 
WRITE (lOT, 7272) 

WRITE (lOT, 1077) 

FORMAT (/'  BRANCH  FROM  NODE  TO  NODE  LENGTH (M)  DIAMETER (CM) ' , 

/) 

DO  1277  J=1,B 

WRITE (lOT, 1255)  J,NF(J) ,NT(J) ,ELL(J) ,D(J) 

CONTINUE 
WRITE (lOT, 7275) 

READ (IN, 7677)  ANS 

IF (ANS  .EQ.  'Y'  .OR.  ANS  .EQ.'y')  GO  TO  111 

IF (ANS  .EQ.  'N'  .OR.  ANS  .EQ.  'n')  GO  TO  1280 

GO  TO  1276 

CONTINUE 

WRITE (lOT, 7281) 

READ  (IN,  1042)  K1 
WRITE (lOT, 7282)  K1 
WRITE (lOT, 1055) 

READdN,  9050)NF(K1)  ,NT(K1)  ,ELL(K1)  ,D(K1) 

WRITE (lOT, 7283) 

READ (IN, 7677)  ANSI 

IF(ANS1  .EQ.  'Y'  .OR.  ANSI  .EQ.  'y' )  GO  TO  19 
IF(ANS1  .EQ.  'N'  .OR.  ANSI  .EQ.  'n')  GO  TO  111 
GO  TO  1285 

DO  28  J=1,B 
ELL(J)=ELL(J)*3.28 
D(J)-D(J)*.3937 
CONTINUE 


DO  33  1=1, B 
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IF  (NF(I)  .LT.  N)  THEN 
A(NF(I)  ,I)=1. 

END  IF 

IF  (NT(I)  .LT.  N)  THEN 
A(NT(I)  ,I)=-1. 

END  IF 
D(I)=D(I)/12 
AREA(I) =PI* (D(I) **2) /4 
LENGTH=ELL(I) 

ADDL=60*D(I) 

ELL ( I ) =LENGTH+ADDL 

33  CONTINUE 

c 
c 
c 
c 

8300  DO  30  K=1,B 

R(K) = (RHO*G*AREA(K) * ( (D(K) ) **2) ) / { 32 . *MU1*ELL (K) ) 
30  CONTINUE 

C 

DO  35  J=1,B 
DO  35  K=1,B 
Y(J,K)  =0. 

35  CONTINUE 


DO  40  1=1, B 

Y(I,I)=R(I) 

40  CONTINUE 

C 
c 

c  QS:  MATRIX  OF  FLOW  SOURCES 

c 

c  S:  NUMBER  OF  FLOW  SOURCES 

c  QB:  SOURCE  BRANCH 

c  QSl:  SOURCE  STRENGTH 

c 

1015  FORMAT  (/'  Input  the  number  of  flow  sources 3X, \) 

1016  FORMAT  (//'  Input  the  branch  of  source' , 13 , 3X, , \) 

1017  FORMAT  (//'  Input  the  strength  of  the  source  (gal/min) : ' , 3X, \) 

9017  FORMAT  (//'  Input  the  strength  of  the  source  (liter/min) : ' , 3X, /) 

WRITE (lOT, 1015) 

READ(IN,1004)  S 
DO  50  1=1, S 

WRITE (lOT, 1016)  I 
READ(IN,1004)  QB 

IFdUNITS  .EQ.  1)  THEN 
WRITE (lOT, 1017) 

READ(IN,1001)  QSl 
ELSE 

WRITE (lOT, 9017) 

READ (IN, 1001)  QSl 
QS1=QS1* .2642 
END  IF 
QSl=QSl/444 
DO  60  K«1,B 
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IF  (K  .EQ.  QB)  THEN 

QS (K, 1) =QS1 

ELSE 

QS(K,1)=0. 

END  IF 

60  CONTINUE 

50  CONTINUE 
c 
c 

C  WRITEdOT, 3012)  (QS(I,1)  ,1  =  1, B) 

C 

c 

C  MATRIX  MANIPULATION  TO  DETERMINE  INDIVIDUAL  BRANCH  FLOW  AND 
c  PRESSURE  DISTRIBUTIONS 
C 

C  AT:  TRANSPOSE  OF  MATRIX  A 

C  AY:  PRODUCT  OF  MATRIX  A  AND  MATRIX  Y 

C  YN:  PRODUCT  OF  MATRIX  AY  AND  MATRIX  AT  (YN=AYAT) 

C 

c 

c 

DO  80  1=1, N1 
DO  80  J=1,B 
80  AT(J, I) =A(I, J) 

C 

C 

C 

CALL  MATMUL(AY,A,Y,N1,B,B) 
c 

CALL  MATMUL(YN,AY,AT,N1,B,N1) 
c 
c 

C  IS:  NODE  FLOW  SOURCE  VECTOR  IS=  -AQS 
C 

CALL  MATMUL(IS,A,QS,N1,B, 1) 

C 

C 

C 

DO  90  1=1, N1 

IS(I,1)=-IS(I,1) 

ISI(I,1)=IS(I,1) 

90  CONTINUE 

C 
C 

c 

C  YNI:  MATRIX  EQUAL  TO  MATRIX  YN  BEFORE  CHOLESKI  INVERSION.  AFTER 
C  DECOMPOSITION  IT  HOLDS  THE  UPPER  TRIANGULAR  MATRIX. 

C  ISI:  MATRIX  EQUAL  TO  MATRIX  IS  BEFORE  INVERSION.  AFTER  INVERSION, 
C  IT  HOLDS  THE  SOLUTION  TO  E=YN(-1)IS 

C 

c 

DO  95  1=1, N1 
DO  95  J=1,N1 

YNI (I, J) =YN(I, J) 

95  CONTINUE 


c 
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CALL  CHOLESKY(YNI, ISI,N1,N1) 
c 
c 

c  E;  NODE  VOLTAGE  MATRIX  E=YN(-1)IS 
C  P:  BRANCH  PRESSURE  P=ATE 

C  Q;  BRANCH  FLOW  RATE  Q=YP+QS 

c  YP:  PRODUCT  OF  MATRIX  Y  AND  MATRIX  P  (YP) 
c 

DO  96  I«1,N1 

E(I,1)=ISI{I,1) 

96  CONTINUE 


CALL  MATMUL(P,AT,E,B,N1,1) 


CALL  MATMUL(YP,Y,P,B,B,1) 

DO  120  1=1, B 

Q(I,1)=YP(I,1)+QS(I,1) 

120  CONTINUE 
CALL  CLS 


REYNOLDS  NUMBER  CALCULATIONS 
DO  5000  1=1, B 

V(I,1)=ABS(Q(I,1)/{AREA(I)  )  ) 

RE (1,1) =rho*V(I, 1) *d{I) /mul 
IF(RE(I,1)  .GT.  2100)  THEN 
WRITE (lOT, 5001)  I 
WRITE (lOT, 5002) 

5001  FORMAT{/'  REYNOLDS  NUMBER  IN  BRANCH' , 13 , 2x, ' EXCEEDS' , \) 

5002  FORMAT {/'  2100.  THE  FLOW  IS  NOT  LAMINAR. ', 2X, \) 

ELSE 

END  IF 

5000  CONTINUE 


1500  FORMAT (3X, 13 , 7X, Fll . 8 , 2X, 2F15 . 5 , IX, FIO . 0) 

1700  FORMAT (3X, 13 , 5X, Fll . 8 , 3X, 2F15 . 5 , 2X, F5 . 0 ) 

1501  FORMAT(///'  Branch' , 4X, ' Flow  Rate  (ft/sec) ', 4x, ' P  in  (psf)',6x, 

+  'P  out  (psf ) ' , 6x, 'Re' , lx, //) 

1601  FORMAT  (///'  Branch' ,  4x,  '  Flow  Rate  (m/sec) ',  5x,  '  P  in  (N/m‘^2)  '  ,  7x, 
+  'Pout  (N/m^2) ' ,6x, 'Re' ,2x,//) 

IFdUNITS  .EQ.  1)  THEN 
WRITE (lOT, 1501) 

DO  7000  1=1, B 

WRITE ( lOT, 1500) I, Q( 1,1) , E { (NF (I) ) , 1) , E { (NT ( I) ) , 1) , RE ( I , 1) 

7000  CONTINUE 

99  FORMAT (/'  WOULD  YOU  LIKE  A  PRINTOUT  OF  THIS  TABLE  (Y/N)?'\) 

98  WRITE (lOT, 99) 

READ (IN, 7677) ANS2 

IF(ANS2  .EQ.  'Y'  .OR.  ANS2  .EQ.  'y')  GO  TO  97 
IF(ANS2  .EQ.  'N'  .OR.  ANS2  .EQ.  'n')  GO  TO  83 
GO  TO  98 
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97  WRITE (IPR, 1501) 

DO  93  1=1, B 

WRITE (IPR, 1500) I,Q(I,1) ,E( (NF(I) ) ,1) ,E( (NT(I) ) ,1) ,RE(I,1) 
93  CONTINUE 

83  CONTINUE 

ELSE 

DO  7470  J=1,N1 
E(J,1)=E{J,1)*47.88 
7470  CONTINUE 

WRITE (lOT, 1601) 

DO  7001  1=1, B 
Q(I,1)=Q(I,1)*.3048 

WRITE (lOT, 1500)  I , Q (I , 1 ) . E ( (NF (I) ) , 1) , E ( (NT ( I ) ) , 1) , RE ( I , 1) 
7001  CONTINUE 

48  WRITE (lOT, 99) 

READ (IN, 7677) ANS2 

IF(ANS2  .EQ.  'Y'  .OR.  ANS2  .EQ.  'y')  GO  TO  47 
IF(ANS2  .EQ.  'N'  .OR.  ANS2  .EQ.  'n')  GO  TO  82 
GO  TO  48 

47  WRITE (IPR, 1601) 

DO  43  1=1, B 

WRITE (IPR, 1500) I, Q( 1,1) , E ( (NF (I) ) , 1) , E ( (NT ( I ) ) , 1) , RE ( I , 1) 
43  CONTINUE 

82  CONTINUE 

END  IF 
END 
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SUBROUTINE  CHOLESKY (A, B, N,NX) 
c 
c 

C  SOLUTION  OF  LINEAR  SYSTEMS  OF  EQUATIONS  BY  THE  CHOLESKI 
C  METHOD  FOR  SYMMETRIC  POSITIVE  DEFINITIVE  MATRICES, 
c 

C  A:  ARRAY  CONTAINING  THE  SYSTEM  MATRIX  (AX=B) 

C  B:  ARRAY  CONTAINING  THE  VECTOR  IF  INDEPENDENT  COEFFICIENTS 

c  C:  AUXILIARY  VECTOR 

C  N:  ORDER  OF  A 

C  NX:  ROW  AND  COLUMN  DIMENSION  OF  A 
c 

DOUBLE  PRECISION  B(NX,1) 

DIMENSION  A (40, 40) 
c 

C  COMPUTE  UPPER  TRIANGULAR  MATRIX  FROM  A  AND  STORE  ALSO  IN  A 
C 

CALL  DECOMP (A, N, NX) 

C 

c  COMPUTE  THE  C  VECTOR  AND  STORE  IN  ARRAY  B 


c 


B(1,1)=B(1,1)/A(1,1) 

c 

WRITE (lOT, 100) 

A(l,l) 

100 

FORMAT (F18. 15) 
DO  10  1=2, N 

c 

WRITE (lOT, 101) 

I,N 

101 

FORMAT (213) 

c 

WRITE (lOT, 102) 

A(I,I)  ,A(N,N) 

102 

FORMAT (2F13 .11) 

D=B(I,1) 

11=1-1 

DO  5  L=1,I1 

5 

D=D-A{L, I) *B(L,1) 

10 

B(I,1)=D/A(I,I) 

B(N,1)=B(N,1)/A(N,N) 

c  COMPUTE  THE  SYSTEM  UNKNOWNS 

AND  STORE  IN  ARRAY  B 

N1=N-1 

DO  30  L=1,N1 

K=N-L 

C 

WRITE (lOT, 105) 

A(K,K)  ,K 

105 

FORMAT (F13.il, 

Kl-K+1 

DO  20  J=K1,N 

13) 

20 

B(K,1)=B(K,1)  -A(K,  J)  *B(J,1) 

30 

B(K,1)=B(K,1)/A(K,K) 

RETURN 

END 
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SUBROUTINE  MATMUL (C, A, B, N, M, L) 
c  This  subroutine  computes  the  matrix  operation  C  =  A  *  B 
c 

c  N:  NUMBER  OF  ROWS  IN  MATRIX  A  AND  C 

c  M:  NUMBER  OF  COLUMNS  IN  A  AND  ROWS  IN  B 

c  L;  NUMBER  OF  COLUMNS  IN  B  AND  C 

c 

DIMENSION  A(40,40)  ,B(40,40)  ,C(40,40) 

DO  20  1=1, N 
DO  20  J=1,L 
C(I,  J)=0.0 
DO  20  K=1,M 

20  C(I,  J)=C(I,  J)+A(I,K)*B(K,  J) 

RETURN 

END 

c 
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c 

SUBROUTINE  DECOMP (A, N, NX) 
c 

C  THIS  SUBROUTINE  PERFORMS  THE  DECOMPOSITION  OF  A  POSITIVE, 

C  DEFINITE,  SYMMETRIC  MATRIX  INTO  AN  UPPER  TRIANGULAR  MATRIX, 
c 

C  A:  ARRAY  ORIGINALLY  CONTAINING  MATRIX  TO  BE  DECOMPOSED. 

C  AT  THE  END  IT  CONTAINS  THE  UPPER  TRIANGULAR  MATRIX. 

C  N:  ORDER  OF  A 

C  NX:  ROW  AND  COLUMN  DIMENSION  OF  A 

INTEGER  N 
DIMENSION  A(40,40) 

IF  (A(l,l)  )1,1,3 

1  WRITE(IOT,2) 

2  FORMAT  ('  ZERO  OR  NEGATIVE  RADICAND' ) 

C  WRITE (lOT, 500)  A(l,l) 

500  FORMAT  ('  A (1 , 1) ' , F13 . 11) 

GO  TO  200 

3  A(1,1)=SQRT(A(1,1)  ) 

DO  10  J=2,N 

10  A(l,  J)=A(1,J)/A(1,1) 

DO  40  1=2, N 

C  WRITE (lOT, 600)  I 

600  FORMAT (13) 

11=1-1 

D=A(I,I) 

DO  20  L=1,I1 

20  D=D-A(L, I) *A(L, I) 

C  WRITE (lOT, 503)  I,D,A(I,I) 

503  FORMAT ( 13 , 2F13 .11) 

IF  (A(I,I) )1,1,21 

21  A{I,I)=SQRT(D) 

IF (I  .EQ.  N)  THEN 
GO  TO  45 

END  IF 

12=1+1 

DO  40  J=I2,N 
D=A(I, J) 

DO  30  L=1,I1 
30  D=D-A(L,I) *A(L, J) 

C  WRITE (lOT, 503)  I,D,A(I,I) 

40  A(I,  J)=D/A(I,  I) 

45  DO  50  1=2, N 
11=1-1 

DO  50  J=1,I1 
50  A(I,J)=0 

200  RETURN 
END 
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APPENDIX  H 


Case  Run  One 


Branch 

Flow  Rate 

Pin 

Pout 

Re 

(ft/sec) 

(psf ) 

(psf) 

1 

.00035966 

- .03996 

.05160 

1299 

2 

.00019242 

.05160 

.01905 

695 

3 

.00015641 

.01905 

-.00741 

565 

4 

.00019242 

-  .00741 

-  .03996 

695 

5 

.00016724 

- .00555 

-  .  03996 

604 

6 

.00016724 

.05160 

.01719 

604 

7 

.00003601 

.01905 

.01164 

130 

8 

- . 00003601 

- .00741 

.00000 

130 

9 

- .00013443 

- .00555 

.01719 

485 

10 

. 00003281 

.01719 

.01164 

118 

11 

.00006882 

.01164 

. 00000 

249 

12 

.00003281 

.00000 

- .00555 

118 

Case  Run  Two 


Branch 

Flow  Rate 
(tn/sec) 

Pin 

(N/sm) 

Pout 

(N/sm) 

Re 

1 

.00013296 

- .28579 

.37451 

698  . 

2 

.00006895 

.37451 

.13925 

362  . 

3 

.00005561 

.13925 

- . 05052 

292  . 

4 

.  00006895 

-  .05052 

- .28579 

362  . 

5 

.  00006401 

-  .04323 

- .28579 

336  . 

6 

.  00006401 

.37451 

.13195 

336. 

7 

.00001333 

.13925 

. 08872 

70  . 

8 

- .00001333 

- .05052 

.00000 

70  . 

9 

-  .00005134 

-  .04323 

.13195 

270  . 

10 

.00001267 

.13195 

. 08872 

67. 

11 

.00002600 

.08872 

.00000 

137  . 

12 

.00001267 

.00000 

-  .04323 

67. 
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